subroutine find_selection_ParticleIDs(x_in,y_in,z_in,r_in,is_final)
  use log
  use myfiles
  use particles
  use params_and_constants
  use hdf5_wrapper
  use indexx
  
  implicit none
  integer*4               :: itype,ipart,file_handle
  integer*4               :: id
  integer*4, allocatable  :: ind(:),ids(:)
  real*4, intent(in)      :: x_in,y_in,z_in,r_in
  real*4                  :: x,y,z,r2,r_in2
  character(len=16)       :: ctype
  logical                 :: is_final
  integer, dimension(0:5) :: Num

  write(logunit,*)'  -identifying particles within selection sphere...'
  write(logunit,*)'   (x,y,z,r) = ',x_in,y_in,z_in,r_in
  call flush(logunit)

  ! Allocate tmp selection arrays
  call allocate_tmp_selection()

  Num(:)         = 0 
  r_in2          = r_in*r_in

  do itype=0,5
     do ipart=0,Header%NumPart_Total(itype)-1

        selectcase(itype) 
        case(0)
           x  = Type0(ipart)%Coordinates(0) - x_in
           y  = Type0(ipart)%Coordinates(1) - y_in
           z  = Type0(ipart)%Coordinates(2) - z_in
        case(1)
           x = Type1(ipart)%Coordinates(0) - x_in
           y = Type1(ipart)%Coordinates(1) - y_in
           z = Type1(ipart)%Coordinates(2) - z_in
        case(2)
           x = Type2(ipart)%Coordinates(0) - x_in
           y = Type2(ipart)%Coordinates(1) - y_in
           z = Type2(ipart)%Coordinates(2) - z_in
        case(3)
           x = Type3(ipart)%Coordinates(0) - x_in
           y = Type3(ipart)%Coordinates(1) - y_in
           z = Type3(ipart)%Coordinates(2) - z_in
        case(4)
           x = Type4(ipart)%Coordinates(0) - x_in
           y = Type4(ipart)%Coordinates(1) - y_in
           z = Type4(ipart)%Coordinates(2) - z_in
        case(5)
           x = Type5(ipart)%Coordinates(0) - x_in
           y = Type5(ipart)%Coordinates(1) - y_in
           z = Type5(ipart)%Coordinates(2) - z_in
        case default
           call abort("invalid particle type")
        end select
        r2 = x*x + y*y + z*z

        if(r2 <= r_in2)then

           selectcase(itype) 
           case(0)
              Type0_tmp(Num(0))%Index            = ipart
              Type0_tmp(Num(0))%ParticleIDs      = Type0(ipart)%ParticleIDs
              Type0_tmp(Num(0))%Coordinates(0:2) = Type0(ipart)%Coordinates(0:2)
           case(1)
              Type1_tmp(Num(1))%Index            = ipart
              Type1_tmp(Num(1))%ParticleIDs      = Type1(ipart)%ParticleIDs
              Type1_tmp(Num(1))%Coordinates(0:2) = Type1(ipart)%Coordinates(0:2)
           case(2)
              Type2_tmp(Num(2))%Index            = ipart
              Type2_tmp(Num(2))%ParticleIDs      = Type2(ipart)%ParticleIDs
              Type2_tmp(Num(2))%Coordinates(0:2) = Type2(ipart)%Coordinates(0:2)
           case(3)
              Type3_tmp(Num(3))%Index            = ipart
              Type3_tmp(Num(3))%ParticleIDs      = Type3(ipart)%ParticleIDs
              Type3_tmp(Num(3))%Coordinates(0:2) = Type3(ipart)%Coordinates(0:2)
           case(4)
              Type4_tmp(Num(4))%Index            = ipart
              Type4_tmp(Num(4))%ParticleIDs      = Type4(ipart)%ParticleIDs
              Type4_tmp(Num(4))%Coordinates(0:2) = Type4(ipart)%Coordinates(0:2)
           case(5)
              Type5_tmp(Num(5))%Index            = ipart
              Type5_tmp(Num(5))%ParticleIDs      = Type5(ipart)%ParticleIDs
              Type5_tmp(Num(5))%Coordinates(0:2) = Type5(ipart)%Coordinates(0:2)
           case default
              call abort("invalid particle type")
           end select
           Num(itype) = Num(itype) + 1

        endif
        
     enddo
  enddo

  write(logunit,*)'   -selection sphere contains:'
  do itype=0,5 
     write(ctype,*)itype
     write(logunit,*)'    -itype = '//trim(adjustl(ctype))//' NumPart = ',Num(itype)
  enddo

  ! Allocate selection arrays
  call allocate_selection()

  ! Sort the selected particles by ParticleIDs
  write(logunit,*)'   -sorting selection particles by index...'
  do itype=0,5

     write(ctype,*)itype

     if(Num(itype) > 0)then 

        write(logunit,*)'    -type = '//trim(adjustl(ctype))
        if(allocated(ind)) deallocate(ind)
        if(allocated(ids)) deallocate(ids)
        allocate(ind(0:Num(itype)-1))
        allocate(ids(0:Num(itype)-1))  

        selectcase(itype)
        case(0)
           ids = Type0_tmp(0:Num(0)-1)%ParticleIDs
           call indexxx(Num(0),ids,ind)
           ind(:) = ind(:) - 1 
           do ipart=0,Num(0)-1
              Type0_Selection(ipart)%ParticleIDs    = Type0_tmp(ind(ipart))%ParticleIDs
              Type0_Selection(ipart)%Index          = Type0_tmp(ind(ipart))%Index
              Type0_Selection(ipart)%Coordinates(:) = Type0_tmp(ind(ipart))%Coordinates(:)
           enddo
        case(1)
           ids = Type1_tmp(0:Num(1)-1)%ParticleIDs
           call indexxx(Num(1),ids,ind)
           ind(:) = ind(:) - 1 
           do ipart=0,Num(1)-1
              Type1_Selection(ipart)%ParticleIDs    = Type1_tmp(ind(ipart))%ParticleIDs
              Type1_Selection(ipart)%Index          = Type1_tmp(ind(ipart))%Index
              Type1_Selection(ipart)%Coordinates(:) = Type1_tmp(ind(ipart))%Coordinates(:)
           enddo
        case(2)
           ids = Type2_tmp(0:Num(2)-1)%ParticleIDs
           call indexxx(Num(2),ids,ind)
           ind(:) = ind(:) - 1 
           do ipart=0,Num(2)-1
              Type2_Selection(ipart)%ParticleIDs    = Type2_tmp(ind(ipart))%ParticleIDs
              Type2_Selection(ipart)%Index          = Type2_tmp(ind(ipart))%Index
              Type2_Selection(ipart)%Coordinates(:) = Type2_tmp(ind(ipart))%Coordinates(:)
           enddo
        case(3)
           ids = Type3_tmp(0:Num(3)-1)%ParticleIDs
           call indexxx(Num(3),ids,ind)
           ind(:) = ind(:) - 1 
           do ipart=0,Num(3)-1
              Type3_Selection(ipart)%ParticleIDs    = Type3_tmp(ind(ipart))%ParticleIDs
              Type3_Selection(ipart)%Index          = Type3_tmp(ind(ipart))%Index
              Type3_Selection(ipart)%Coordinates(:) = Type3_tmp(ind(ipart))%Coordinates(:)
           enddo
        case(4)
           ids = Type4_tmp(0:Num(4)-1)%ParticleIDs
           call indexxx(Num(4),ids,ind)
           ind(:) = ind(:) - 1 
           do ipart=0,Num(4)-1
              Type4_Selection(ipart)%ParticleIDs    = Type4_tmp(ind(ipart))%ParticleIDs
              Type4_Selection(ipart)%Index          = Type4_tmp(ind(ipart))%Index
              Type4_Selection(ipart)%Coordinates(:) = Type4_tmp(ind(ipart))%Coordinates(:)
           enddo
        case(5)
           ids = Type5_tmp(0:Num(5)-1)%ParticleIDs
           call indexxx(Num(5),ids,ind)
           ind(:) = ind(:) - 1 
           do ipart=0,Num(5)-1
              Type5_Selection(ipart)%ParticleIDs    = Type5_tmp(ind(ipart))%ParticleIDs
              Type5_Selection(ipart)%Index          = Type5_tmp(ind(ipart))%Index
              Type5_Selection(ipart)%Coordinates(:) = Type5_tmp(ind(ipart))%Coordinates(:)
           enddo
        endselect
        deallocate(ind)
        deallocate(ids)
     endif
  enddo

  ! Free up some memory
  call deallocate_tmp_selection()

 if(is_final)then
    NumSelected = Num
 else
    NumGathered = Num
 endif

  write(logunit,*)'  -done'
  call flush(logunit)

end subroutine find_selection_ParticleIDs

! **********************************************************

subroutine trace_selection_indices()
  use log
  use myfiles
  use particles
  use params_and_constants
  use indexx
  use seek_module

  implicit none
  integer*4              :: itype,ipart,file_handle
  integer*4              :: id,this_ind
  integer*4, allocatable :: ind(:),ids(:)
  real*4                 :: x,y,z,r2
  character(len=16)      :: ctype

  write(logunit,*)'  -finding indices of selected particles...'
  call flush(logunit)

  ! Allocate global for traced particles
  call allocate_traced()

  ! Sort each particle type separately 
  NumTraced(:) = 0
  do itype=0,5

     if(NumSelected(itype) == 0) cycle

     write(ctype,*)itype     
     write(logunit,*)'   -type = '//trim(adjustl(ctype))

     ! Allocate local array for sorting particle ids from this snapshot
     if(allocated(ind)) deallocate(ind)
     allocate(ind(0:Header%NumPart_Total(itype)-1))

     ! Also allocate local array for accepting the ParticleIDs - remember that
     ! Nehalems (e.g. shrek301.swin.edu.au) don't like passing structure elements
     if(allocated(ids)) deallocate(ids)
     allocate(ids(0:Header%NumPart_Total(itype)-1))     

     ! Sort by ParticleIDs
     write(ctype,*)itype     
     write(logunit,*)'    -sorting snapshot particles by index'

     selectcase(itype)
     case(0)
        ids = Type0%ParticleIDs
     case(1)
        ids = Type1%ParticleIDs
     case(2)
        ids = Type2%ParticleIDs
     case(3)
        ids = Type3%ParticleIDs
     case(4)
        ids = Type4%ParticleIDs
     case(5)
        ids = Type5%ParticleIDs
     endselect

     call indexxx(Header%NumPart_Total(itype),ids,ind)

     selectcase(itype)
     case(0)
        Type0(:)%Index = ind(:) - 1
     case(1)
        Type1(:)%Index = ind(:) - 1
     case(2)
        Type2(:)%Index = ind(:) - 1
     case(3)
        Type3(:)%Index = ind(:) - 1
     case(4)
        Type4(:)%Index = ind(:) - 1
     case(5)
        Type5(:)%Index = ind(:) - 1
     endselect
     
     ! Now we can cross-reference with recursive bisection. Some things to note:
     ! 1 - all gas and cdm (incl. boundary) particles should be found, as they either
     !     disappear or are conserved as we go forwards in time.
     ! 2 - stars *appear* over time, so disappear when tracing back. Hence we must look
     !     for any unfound stars in the gas arrays at this snapshot.
     write(logunit,*)'    -cross-referencing selected particles with snapshot'
     do ipart=0,NumSelected(itype)-1
        
        selectcase(itype)
        case(0)
           this_ind = rb_seek(0,Type0_Selection(ipart)%ParticleIDs)
           if(this_ind == -1) call abort(" could not trace Type0 particle ")
           Type0_Traced_Index(NumTraced(0)) = Type0(this_ind)%Index
           NumTraced(0) = NumTraced(0) + 1
        case(1)
           this_ind = rb_seek(1,Type1_Selection(ipart)%ParticleIDs)
           if(this_ind == -1) call abort(" could not trace Type1 particle ")
           Type1_Traced_Index(NumTraced(1)) = Type1(this_ind)%Index
           NumTraced(1) = NumTraced(1) + 1
        case(2)
           this_ind = rb_seek(2,Type2_Selection(ipart)%ParticleIDs)
           if(this_ind == -1) call abort(" could not trace Type2 particle ")
           Type2_Traced_Index(NumTraced(2)) = Type2(this_ind)%Index
           NumTraced(2) = NumTraced(2) + 1
        case(3)
           this_ind = rb_seek(3,Type3_Selection(ipart)%ParticleIDs)
           if(this_ind == -1) call abort(" could not trace Type3 particle ")
           Type3_Traced_Index(NumTraced(3)) = Type3(this_ind)%Index
           NumTraced(3) = NumTraced(3) + 1
        case(4)

           ! Check if even worth looking in the stars first 
           if(Header%NumPart_Total(4) > 0)then 

              this_ind = rb_seek(4,Type4_Selection(ipart)%ParticleIDs)
              if(this_ind == -1) then ! Don't abort yet, look in the gas!!!

                 this_ind = rb_seek(0,Type4_Selection(ipart)%ParticleIDs)
                 if(this_ind == -1) call abort(" could not trace Type4 particle (even in the gas) ")

                 ! So we've found it in the gas, therefore attribute as such 
                 Type0_Traced_Index(NumTraced(0)) = Type0(this_ind)%Index
                 if(ipart == 0)write(logunit,*)ipart,Type0_Traced_Index(NumTraced(0)),&
                      Type4_Selection(ipart)%ParticleIDs,Type0(Type0_Traced_Index(NumTraced(0)))%ParticleIDs
                 NumTraced(0) = NumTraced(0) + 1
              else                 

                 ! Found in stars
                 Type4_Traced_Index(NumTraced(4)) = Type4(this_ind)%Index
                 if(ipart == 0)write(logunit,*)ipart,Type4_Traced_Index(NumTraced(4)),&
                      Type4_Selection(ipart)%ParticleIDs,Type4(Type4_Traced_Index(NumTraced(4)))%ParticleIDs
                 NumTraced(4) = NumTraced(4) + 1
              endif

           else

              ! Just look in the gas
              this_ind = rb_seek(0,Type4_Selection(ipart)%ParticleIDs)
              if(this_ind == -1) call abort(" could not trace Type4 particle (even in the gas) ")
              ! So we've found it in the gas, therefore attribute as such 
              Type0_Traced_Index(NumTraced(0)) = Type0(this_ind)%Index
              NumTraced(0) = NumTraced(0) + 1

           endif

        case(5)
           this_ind = Type5(rb_seek(itype,Type5_Selection(ipart)%ParticleIDs))%Index
           if(this_ind == -1) call abort(" could not trace Type5 particle ")
           Type5_Traced_Index(NumTraced(5)) = Type5(this_ind)%Index
           NumTraced(5) = NumTraced(5) + 1
        endselect
        
     enddo

  enddo

  ! Free up memory 
  deallocate(ind)
  deallocate(ids)

  ! Report evolution of traced version of the selection particle numbers
  write(logunit,*)'   -particle selection traced back to:'
  write(logunit,*)'    -type0: selection = ',NumSelected(0),' traced = ',NumTraced(0)
  write(logunit,*)'    -type1: selection = ',NumSelected(1),' traced = ',NumTraced(1)
  write(logunit,*)'    -type2: selection = ',NumSelected(2),' traced = ',NumTraced(2)
  write(logunit,*)'    -type3: selection = ',NumSelected(3),' traced = ',NumTraced(3)
  write(logunit,*)'    -type4: selection = ',NumSelected(4),' traced = ',NumTraced(4)
  write(logunit,*)'    -type5: selection = ',NumSelected(5),' traced = ',NumTraced(5)
  write(logunit,*)'   -done'
  call flush(logunit)

end subroutine trace_selection_indices

! **********************************************************

subroutine compute_traced_com()
  use log
  use myfiles
  use particles
  use params_and_constants
  use indexx
  use seek_module

  implicit none
  integer*4              :: itype,ipart
  integer*4              :: ind
  real*4                 :: MassTot

  write(logunit,*)'  -finding centre of mass of traced particles...'
  call flush(logunit)

  ! Sort each particle type separately 
  COM(:)  = 0
  MassTot = 0.0e0
  do itype=0,5
     if(NumTraced(itype) == 0) cycle
     do ipart=0,NumTraced(itype)-1
        selectcase(itype)
        case(0)     
           ind     = Type0_Traced_Index(ipart)
           COM(0)  = COM(0)  + Type0(ind)%Coordinates(0)*Type0(ind)%Mass           
           COM(1)  = COM(1)  + Type0(ind)%Coordinates(1)*Type0(ind)%Mass
           COM(2)  = COM(2)  + Type0(ind)%Coordinates(2)*Type0(ind)%Mass
           MassTot = MassTot + Type0(ind)%Mass
        case(1)     
           ind     = Type1_Traced_Index(ipart)
           COM(0)  = COM(0)  + Type1(ind)%Coordinates(0)*Type1(ind)%Mass           
           COM(1)  = COM(1)  + Type1(ind)%Coordinates(1)*Type1(ind)%Mass
           COM(2)  = COM(2)  + Type1(ind)%Coordinates(2)*Type1(ind)%Mass
           MassTot = MassTot + Type1(ind)%Mass                 
        case(2)                               
           ind     = Type2_Traced_Index(ipart)                 
           COM(0)  = COM(0)  + Type2(ind)%Coordinates(0)*Type2(ind)%Mass           
           COM(1)  = COM(1)  + Type2(ind)%Coordinates(1)*Type2(ind)%Mass
           COM(2)  = COM(2)  + Type2(ind)%Coordinates(2)*Type2(ind)%Mass
           MassTot = MassTot + Type2(ind)%Mass                 
        case(3)                               
           ind     = Type3_Traced_Index(ipart)                 
           COM(0)  = COM(0)  + Type3(ind)%Coordinates(0)*Type3(ind)%Mass           
           COM(1)  = COM(1)  + Type3(ind)%Coordinates(1)*Type3(ind)%Mass
           COM(2)  = COM(2)  + Type3(ind)%Coordinates(2)*Type3(ind)%Mass
           MassTot = MassTot + Type3(ind)%Mass                 
        case(4)                               
           ind     = Type4_Traced_Index(ipart)                 
           COM(0)  = COM(0)  + Type4(ind)%Coordinates(0)*Type4(ind)%Mass           
           COM(1)  = COM(1)  + Type4(ind)%Coordinates(1)*Type4(ind)%Mass
           COM(2)  = COM(2)  + Type4(ind)%Coordinates(2)*Type4(ind)%Mass
           MassTot = MassTot + Type4(ind)%Mass                 
        case(5)                                
           ind     = Type5_Traced_Index(ipart)                
           COM(0)  = COM(0)  + Type5(ind)%Coordinates(0)*Type5(ind)%Mass           
           COM(1)  = COM(1)  + Type5(ind)%Coordinates(1)*Type5(ind)%Mass
           COM(2)  = COM(2)  + Type5(ind)%Coordinates(2)*Type5(ind)%Mass
           MassTot = MassTot + Type5(ind)%Mass
        endselect
     enddo
  enddo

  write(logunit,*)'  -MassTot = ',MassTot
  COM(:) = COM(:) / MassTot
  write(logunit,*)'  -COM = ',COM

  write(logunit,*)'   -done'
  call flush(logunit)

end subroutine compute_traced_com
